; docformat = 'rst'
;
;+
;
; :Purpose:
;   Calculate CO2-equivalent emissions or mole fractions
;
; :Inputs:
;   input: Either emissions (default) or mole fractions (when /mole_fraction is set)
;
; :Requires:
;   mr_GWP.csv file containing global warming potentials
;   mr_radiative_efficiency.csv if mole fraction CO2e is requested
;
; :Keywords:
;   mole_fraction: set to true (1) if the CO2e mole fracion is required
;   reverse: Undo GWP scaling (i.e. calculate real emissions from CO2e emissions).
;     For emissions only.
;   
; :Outputs:
;   CO2e emissions or CO2e mole fractions
;
; :Example::
;   emissions_CO2e=mr_CO2e(emissions, 'CH4')
;   mole_fraction_CO2e=mr_CO2e(mole_fraction, 'CH4', /mole_fraction)
;
; :History:
;   Written by: Matt Rigby, University of Bristol, Jun 14, 2013
;
;-
function mr_co2e, input, pollutant, mole_fraction=mole_fraction, reverse=reverse, sar=sar
  
  if ~keyword_set(mole_fraction) then begin

    ;Get emissions CO2 equivalent
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    data=read_csv(file_which('mr_GWP.csv'), record_start=2)
    pollutantRef=reform(data.(0))
    if keyword_set(sar) then begin
      GWPref=float(reform(data.(2)))
    endif else begin
      GWPref=float(reform(data.(1)))
    endelse

    wh=where(pollutantRef eq pollutant, count)
    if count gt 0 then begin
      GWP=GWPref[wh]
    endif else begin
      message, "Can't find " + pollutant + " in mr_GWP.csv"
    endelse
  
    if ~keyword_set(reverse) then begin
      output=input*GWP[0]      
    endif else begin
      output=input/GWP[0]      
    endelse

  endif else begin

    ;Get mole fraction CO2 equivalent
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    rf=mr_radiative_forcing(input, pollutant)
    CO2_rf=mr_radiative_forcing(400., 'CO2')
    
    ;Get CO2 RF parameters
    data=read_csv(file_which('mr_radiative_efficiency.csv'), record_start=2)
    pollutantRef=reform(data.(0))
    alphaRef=float(reform(data.(1)))
    C0Ref=float(reform(data.(2)))
    wh=where(pollutantRef eq 'CO2', count)
    if count gt 0 then begin
      alpha=alphaRef[wh]
      C0=C0Ref[wh]
    endif else begin
      message, "Can't find CO2 in mr_radiative_efficiency.csv"
    endelse

    output=replicate(!values.f_nan, (size(input))[1:(size(input))[0]])
    wh=where(finite(input), count)
    if count gt 0 then begin
      output[wh]=C0[0]*exp((rf[wh] + CO2_rf)/1000./alpha[0]) - 400.
    endif
    
  endelse

  return, output

end